Short range repulsive interatomic interactions in energetic processes in solids 
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The repulsive interaction between two atoms at short distances is studied in order to explore 
the range of validity of standard first-principles simulation techniques and improve the available 
short-range potentials for the description of energetic collision cascades in solids. Pseudopotentials 
represent the weakest approximation, given their lack of explicit Pauli repulsion in the core-core 
interactions. The energy (distance) scale realistically accessible is studied by comparison with all- 
electron reference calculations in some binary systems. Reference calculations are performed with 
no approximations related to either core (frozen core, augmentation spheres) or basis set. This is 
important since the validity of such approximations, even in all-electron calculations, rely on the 
small core perturbation usual in low-energy studies. The expected importance of semicore states is 
quantified. We propose a scheme for improving the electronic screening given by pseudopotentials 
for very short distances. The results of this study are applied to the assessment and improvement 
of existing repulsive empirical potentials. 

PACS numbers: 71.15.-m,71.15.Nc,79.20.Ap,61.80.-x,34.20.Cf 
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I. INTRODUCTION 

A high-energy ion (tens to hundreds of keV) propa- 
gating through the system can approach other atoms to 
very close distances, typically between 1 A and 0.1 A. 
The strong repulsive potential determines the scattering 
process. Ion implantation, laser ablation and displace- 
ment cascades produced by a-decay events in materials 
used for nuclear waste immobilization, are examples in 
which these energetic collisions play an important role. 
The energetic ion loses energy due to inelastic collisions 
with the electrons and other atoms in the materiali. The 
latter, called nuclear stopping, dominates for relatively 
low energies (^ 1 keV/amu), while electronic stopping 
dominates for higher energies. 

Nuclear stopping processes represent an enormous 
challenge because the kinetic energies involved allow for 
displacements of the atoms in large regions, and the en- 
ergy dissipation takes place over long time scales (at least 
tens of picoseconds). Complicated structural configu- 
rations are obtained, with atoms breaking and making 
new bonds, and ionisation processes affecting the inter- 
actions between the particles. The complexity in physical 
and chemical interactions, the large sizes, and the long 
times involved in the problem demands high-efliciency 
first-principles calculations. 

Many Density Functional Theory (DFT)^ implementa- 
tions used by the condensed matter community are based 
in the pseudopotential approximation, in which the core 
electrons are considered as frozen in their atomic config- 
uration, and only the valence electrons are responsible of 
the chemical and physical properties of the solid, result- 
ing in a substantial reduction of the degrees of freedom to 
be solved and hence in much faster calculations than the 
equivalent all-electron. For the very challenging problems 
in this area, it would be highly desirable to asses the use- 
fulness and the range of validity of such approximation. 
It has to be remembered, however, that conventional 



all-electron approaches are not much better suited for 
very short distances. If augmented plane waves (APW) 
are used to solve numerically the all-electron problem in 
the core region, problems arise when the augmentation 
spheres overlap. The need for an accurate basis set for 
the core electrons might prevent the use of available gaus- 
sian basis sets, normally devised for unperturbed core 
electrons. 

Reference calculations are thus needed to test the accu- 
racy of pseudopotential-based first principles calculations 
for describing electronic screening at short distances. Ex- 
perimental uncertainties in the determination of the po- 
tentials are usually of the order of 10% or more. We 
have chosen to follow the approach of Nordlund et. aL-, 
and use as a reference calculations at the Hartree-Fock 
limit, where all these problems are eluded and the only 
remaining approximations are the absence of electronic 
correlation and of relativistic effects. In this work we 
have studied the relative importance of each of the ef- 
fects that could be relevant to the problem, and assessed 
the validity of the pseudopotential approximation. 



II. METHODOLOGY 

For high energies the interatomic separation in a head- 
on collision can be small, making the binary collision 
(BC) an appropriate approximation. In this case, the 
complexities of the repulsive interaction between atoms 
can be modelled by a simple two-body interatomic poten- 
tial for short distances. The interatomic potential V{R) 
is defined as the difference between the total energy of 
the two atoms at a distance R, and the energy of the 
isolated atoms: 



V{R) ^ E{R) ~ E{oo) 



(1) 



We use two descriptions to characterize this potential: 
6V{R), and the screening function $(i?). The former is 
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defined as 



6V{R) = V{R) - 



R 



(2) 



where Zi and Z2 are the atomic numbers of the two 
atoms. This can be used to characterize the potentials at 
r — > 0. In the hmit, we have the joint-atom solution, with 
just one atom with atomic number Z\ + Zi. Notice that 
(5F(0) = E''^°-^{Zx + Z2) " - £;"*°"(Z2), and 

we can use all-electron atomic calculations (without the 
problems related to basis commented before) to obtain 
information of this limit of the repulsive potential. 

A different recasting of the interatomic potential that 
is also very useful**, is obtained if we consider the effect of 
screening due to the electronic cloud around the nuclei: 



(3) 



where $(i?) is defined as the screening function. For 
very short distances, there is almost no screening of the 
internuclear repulsion because the probability of finding 
electrons in the interatomic region is very small, and 
$(0) — > 1. For long distances the potential goes to 
zero faster than the Coulomb potential and $(c!o) — > 
(for neutral atoms). Some popular parametrizations 
of this potential are the ones due to Moliere, Lenz- 
Jenson, Born-Mayer, and the one by Ziegler, Biersack 
and Littmark^ (ZBL). We will compare our results with 
the latter, that takes the form: 



(4) 



where the numerical coefficients and hi are fitted to 
a variety of interatomic repulsive potentials for different 
pairs. This gives a relatively good agreement with ex- 
periment, specially for light elements. Nevertheless, the 
deviation with experimental ion scattering^ data can be 
as large as 20%. 

The DFT energy i?(i?) is obtained with the Siesta 
method^i using separable^ norm-conserving TrouUier- 
Martins pseudopotentialsiS. The valence wave func- 
tions are expanded in pseudoatomic numerical orbitals 
including multiple-^ and polarization function. Periodic 
boundary conditions require the definition of a supercell 
large enough for the interaction between replicas of the 
atoms to be sufficiently small (L~20A) . An uniform real- 
space grid defined by an equivalent plane- wave cutoff of 
350 Ry was used for numerical integration. 

In our simulations, scalar relativistic pseudopotentials 
are generated using different electronic configurations, 
including semicore states, partial core corrections, and 
different exchange-correlation functionals. The use of 
semicore shells, that is, core electrons considered as part 
of the valence, has the benefit of increasing the ac- 
curacy as compared to all-electron calculations, mak- 
ing the pseudopotential harder, and more demanding 



computationally. In table 13 the different configurations 
used for core/ valence electrons are shown. Partial core 
correctionsii (PCC) account at least partially for non- 
linearity of the exchange-correlation functional due to 
the charge-density overlap between core and valence elec- 
trons. Notice that including scalar relativistic effects, as 
well as the correlation, represents, at least in some re- 
spects, an improvement to HF. 

We compare the potentials obtained with Siesta with 
those obtained using the same numerical Hartree-Fock 
method as Nordlund et. alA and developed by Laaksonen 
et. alMi. This HF approach, uses wave-functions of the 
general form: 



(27r)i/2 



u(^,77) 



(5) 



and the Hartree-Fock equations are reduced to partial 
differential equations for the function u(^,7y) in prolate 
spheroidal coordinates. The equations are then dis- 
cretized using finite differences. 

This numerical method is used as a reference in our 
electronic structure calculation. The effect of electronic 
correlation was shown not to be relevant, at least for the 
light element potential studied by Nordlund et. al^. We 
confirm this result on other systems. We also checked 
that relativistic effects not included in the Hartree-Fock 
method, are important but not relevant for the purpose 
of this work. Thus, albeit the computational cost of HF 
restricts its application to relatively simple systems, it 
appears to be a good tool to gauge our approximation 
and test the accuracy achieved with other methods. 



III. RESULTS 

The pairs considered in this work were C-C, 0-0, Si- 
Si, Ca-0, and Ca-Ca with both HF and DFT, and U-0 
with only DFT. The screening functions <&(i?) are shown 
in figures m and 13 Deviations from HF become relevant 
for distances R > 0.7A, when the attractive bonding in- 
teraction begins to be important. It is in this region 
of bonding character in the interatomic potential where 
the different descriptions of exchange and correlation give 
different results. It is, however, not important for this 
study. We used the Ceperley-Alder-i? parametrization of 
exchange-correlation in the Local Density Approximation 
(LDA) , and the Perdew, Burke and Ernzerhof— scheme 



TABLE I: Electronic configurations considered for the pseu- 
dopotential generation. 
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FIG. 1: (Color online) Screening function for 0-0 with differ- 
ent approximation for exchange-correlation (spin unpolarized 
LDA, spin polarized LDA and GGA, HF and ZBL). The ref- 
erence configuration for the pseudopotential generation was 
2s^2p'' with core radius rc = O.6IA, except for the hard LDA, 
where rc = 0.53 A. The inset shows the ratio between the 
difference in computed potentials and the HF result. 

for Gradient-corrected functional (GGA), both with and 
without spin polarization. We observed (Fig. 0) that 
any of these approximations for exchange and correla- 
tion gives essentially the same repulsive potential. The 
agreement with HF is good in general, even for short 
distances, showing that electronic correlation plays only 
a minor role in the repulsive potential. Note that ZBL 
gives a good description of the repulsive HF potential for 
light elements (C, O). For larger elements this agreement 
is not so good (see below). 

A. Analysis of the different approximations 

1. Effect of semicore states 

In Fig. |21we analyze the influence of semicore states in 
the screening. When the core electrons of one atom be- 
gin to overlap with the core electrons of the other atom, 
the Pauli repulsion between them is not properly taken 
into account in the pseudopotential approximation and 
consequently the result deviates from the Hartree-Fock 
description. The use of semicore states in the pseudopo- 
tential improves this, with part of the core included as 
valence. This can clearly be seen in the top Fig. 2a, com- 
paring the HF results with those obtained with pseudopo- 
tentials. In 2b, we observe the influence of these semicore 
states in the hardness/softness of the repulsive potential. 
The absence of semicore electrons in the valence, makes 
V{R) softer. Fig. 2c shows the deviation of the potential 
with respect to the HF value. Although this result was 
expected, it is important to quantify it with respect to 



other effects, and to obtain an easy rule to assess the 
limit of validity of the pseudopotential approximation. 

When the interatomic distance is comparable to bond- 
ing distances (R~lA), the valence electrons are spread 
around the atoms and in the bonding region, screening 
the Coulomb interaction between the nuclei. As the dis- 
tance between the nuclei is reduced, the repulsive elec- 
tronic interactions (Coulomb and Pauli) in the bonding 
region increase, and the electrons are more unlikely to be 
in the confined area between the nuclei. This is shown 
in Fig. 121 where the valence charge distribution around 
the atoms is displaced from the interatomic center, to 
the surrounding area as the distance decreases. In a sim- 
ulation within the pseudopotential approximation, the 
valence electrons have a smaller participation in the nu- 
clear screening, but the core electrons are assumed to be 
fixed in their initial configuration, and thus, their con- 
tribution to the screening is roughly the same. For this 
reason $(-R) is always smaller than the expected screen- 
ing function with an all electron approach. 

Core electrons deep in the inner shells (see Fig. 2d), are 
more strongly bonded to the nuclei and their participa- 
tion in the collision would only be relevant for very high 
collision energies when the interatomic distance becomes 
really small. For that energy scale, the electronic stop- 
ping power is more important than the nuclear stopping 
power. Hence, an accurate description of the repulsive 
potential for shorter (higher) distances (energies), is be- 
yond the scope of this study. 

Only the outer core-shell electrons are included in 
the semicore. How localized these electrons are can be 
roughly estimated using the spread of the corresponding 
atomic Kohn-Sham orbitals in a DFT all-electron calcula- 
tion. In tabled] we show for each atomic shell, the radii 
rgo and rgg, in which 90% or 99% of the atomic wave- 
function's norm is localized. Two shells corresponding 
to different atoms will likely interact, if the interatomic 
distance is smaller than sum of the radii at which the 
charge is more localized. According to this, and consid- 
ering rgo as a reference, in a Si-Si collision, we can say 
that 2s electrons will interact with 2s or 2p valence elec- 
trons in the other atom for distances R^O.QSA, which 
is an appreciable distance (see arrows in Fig. 2a). So, 
if we want to have a good description of the collision for 
shorter ranges, we would need to include 2s electrons into 
the semicore. With that, we would ensure a reasonable 
description of the interactions up to distances of the or- 
der of O.2A, when the Is electrons start to play a role. 
This would explain why the pseudopotential without any 
semicore state gives such a soft repulsive potential. 

The same argument may be used for Is electrons in 
a C-C collision for distances smaller than 0.4A, and for 
the 2s-2p electrons in a Ca-Ca collision below ~0.6A. 
In the Ca-Ca case, both the 3s and 3p shells have to be 
included in the semicore for reasonable screening down to 
2rgo, in accordance to the model. In Si-Si, however, the 
2p electrons seem to be enough, and the 2s apparently do 
not affect the screening for distances smaller than ~lA. 
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FIG. 2: (Color online), (a) Screening function for C-C, Si-Si and Ca-Ca, with different electronic configurations for the 
pseudopotential. (b) Interatomic repulsive potential V{R) in keV. (c) Deviation of the potential with respect to the HF limit, 
(d) Charge distribution, r^p{r), showing the core electronic shells. (*)3p®4s^ differs from 3p^4s^ just in the number of ("'s used 
to describe the 3p orbitals (double-^", and single-C^ respectively). The small arrows in (a) show the distance at which electronic 
clouds from the core start to overlap with each other (R=2r9o, see table llT^ . 



2. Effect of the core radius 

In the pseudopotential approximation, it is assumed 
that the core electrons of different atoms do not overlap. 



For the interatomic distances of interest in this work, this 
is not the case. Although the core radii of the pseudopo- 
tentials (re) play an important role in the description of 
ground state properties of real materials, the effects for 
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TABLE II; Radii (in A) containing 90% or 99% of the total charge for each atomic orbital for C, O, Si and Ca. Core electron 
eigenenergies (e in eV) are also shown as computed with an all-electron atomic DFT, as well as experimental values for the 
corresponding binding energiesM . 
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the repulsive potential at short distances is not as dra- 
matic. Changes of up to 40% in Tc give differences in 
the repulsive potentials smaller than 2% when the inter- 
atomic distance is smaller than 2rc, whereas differences 
can be as large as 20% when the distance is close to 2rc. 



3. Effect of charge state, correlation, and relativistic 
corrections 



that we use as a reference. With an all-electron DFT- 
based atomic calculation which includes relativistic ef- 
fects (scalar and spin-orbit), we use SV{0) to estimate 
an upper limit of the relativistic corrections. As an ex- 
ample a limit of ~1.3 keV is obtained in a Ca-Ca col- 
lision. Although this is a considerable correction, the 
energy scales we are interested in (see Fig. 0}, are of the 
order of tens of keV. Furthermore, the relativistic effects 



The charge state of the shooting ion can play an impor- 
tant role in the description of the electronic energy loss 
for intermediate energies. For the short range interatomic 
repulsion, however, these effects are not so relevant: since 
electrons move out of the internuclear region for short in- 
teratomic distances (Fig. O, there will be no important 
effect on the repulsive potential if the ionization involved 
not too deep electrons. 

Correlation and relativistic corrections are known to 
be relevant for very short distances or very heavy nu- 
clei. In the pseudopotential method, scalar relativistic 
corrections are included, but not in the HF approach 
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FIG. 3; Contour plot for the charge distribution at different 
interatomic distances for the Ca-Ca pair (0.1, 0.3, 0.5 and lA, 
going from left to right and up to down). 
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FIG. 4: (Color online) Non-coulombic interaction, SV{R), for 
C-C (top) and Ca-Ca (bottom), and comparison with rela- 
tivistic effects. The solid ticks at ~3.3 keV (for C-C) and 
~60 keV (Ca-Ca) represent the atomic energy of Mg and Zr 
(atomic numbers twice as big as C and Ca) with an all-electron 
atomic DFT. The dashed tickmark represents the atomic en- 
ergy with relativistic corrections (for C if falls over the solid 
lines in the energy scale used). 



6 



0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 




I I I I I I I I I I I I 

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 

Interatomic distance R (A) 



FIG. 5: (Color online) Screening function, ^{R), and repul- 
sive interatomic potential, V{R) in the inset, for the U-O 
pair. The arrow shows the radius at which the n=5 shell of 
U overlaps with the n=l shell of O. Bottom figure represents 
the electronic charge distribution for U, as computed with an 
atomic all electron DFT method. 

are expected to be smaller for distances ^0.4A. The dif- 
ference between DFT and HF (both relativistic or both 
non-relativistic) is of the order of ~5 eV and gives an 
idea of the magnitude of correlation effects. It is evident 
from the figure that the approximations involved in the 
pseudopotential approach are much more relevant than 
correlation or relativistic effects. 



B. Practical realization: Heavy ions 

We can now try to obtain some qualitative informa- 
tion about the interatomic potential for heavy ions. The 
use of pseudopotentials for equilibrium (low energy) sit- 
uations of lanthanide and actinide compounds has been 
shown to be accurateiS^ii when the effects coming from 
overlap between core and valence are properly treated 
(either with the inclusion of semicore states, or with par- 
tial core correction), even if only scalar relativistic cor- 
rections are considered. Relativistic spin-orbit terms can 
also be important in these elements, and the results ob- 
tained within this method should be taken with caution. 

We have computed the uranium-oxygen interatomic 
potential V{R) with different pseudopotentials, and com- 
pared it to ZBL (see Fig. [SJ. The electronic configura- 
tion is shown in table and the localization radii rgg 
and rgg , in table IIIII Using the same arguments as be- 



fore, we could expect that using a semicore with 6s and 
6p electrons in the valence, the U-O interatomic poten- 
tial obtained would be appropriate up to distances of the 
order of O.TA. Note that the ZBL potential would give 
a lower screening than expected (at least for distances 
larger than 0.7 A), giving a harder description of the bi- 
nary collision. 

C. Correction of the screening function 

A systematic feature of the pseudopotential calcula- 
tions is that in the limit of small distances, the screen- 
ing functions obtained do not approach unity. This is 
due to the frozen core approximation and the fact that 
the inner shells of the electronic cloud are not contribut- 
ing to the screening. We propose a method improve 
by adding a continuous function that gives the proper 
limit ($(0) — > 1) for values below the radius at 
which the core states begin to overlap. By inspection 
of figures |21 and |S1 we could use the simplest function 
$(i?) = e"^ with parameters a and /3 fixed by the 
conditions over <I>(ro) and ^{ro), giving continuity of the 
screening and its derivative (and thus, the interatomic 
force) at r-g. Notice that the matching radii would de- 
pend on the particular configuration used for the genera- 
tion of the pseudopotential. We propose to use the sums 
of the rgg for the relevant core states of both atoms con- 
sidered: rg = rgQ°™^ -|-rgQ°™^. In figurelHIwe compare the 
previous screening functions with the ones corrected with 
this scheme. This will allow realistic molecular dynamic 
simulations with pseudopotentials, even if the atoms get 
extremely close in a high energy collision. Numerical fits 
of these potentials have already been used to do semiem- 
pirical simulations of collision cascades for studies of ra- 
diation damage in CaTiOgiSi. 

IV. CONCLUSIONS 

The validity of the pseudopotential approximation has 
been checked at close distances as found in high ener- 
getic collisions in materials. We have compared the inter- 
atomic repulsive potential computed with an all-electron 
Hartree-Fock method, with that obtained with pseudopo- 
tentials. In this HF reference method, no approximation 
is assumed with respect to the core shells, or the basis 
used to represent the electronic states. The validity of 
the HF approximation itself has been assessed in order 
to use it as a reference for the pseudopotentials. The 
effect of different approaches to treat the electronic cor- 
relation was previously shown"^ not to give substantial 
differences in the screening function, a result confirmed 
in this study. The relativistic corrections, though im- 
portant, can be considered small for the energy scales of 
interest in this work. 

It has been shown that as the nuclear distance de- 
creases, the electrons are expelled from the internuclear 
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region to the surrounding area, and the screening is com- 
pletely dominated by the deeper electrons, bound at 
short length scales. If the core is frozen at length scales 
longer or comparable to the internuclear distance, the 
screening is larger than expected and the pseudopoten- 
tial description fails. We have shown this effect to be the 
most important one in the problem, actually the only one 
of relevance at the energy scales involved. 

We have quantified to what extent the use of deeper 
electrons as part of the valence configuration is essential 
to describe the repulsive potential for relatively short in- 
ternuclear distances. The localization of the core elec- 
trons, rgo, has been found as the most appropriate cri- 
terion used to decide which states have to be included 
for a particular range of collision energies. The pseu- 
dopotential description would then be valid until the 
more external core shells not included in the calculation 



start to overlap. Finally, we have proposed a scheme 
to correct the screening given in the pseudopotential ap- 
proximation, that recovers the right trend for very short 
distances, and improves the use of pseudopotentials for 
the description of energetic collisions with first principles 
simulations. 
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TABLE III: Estimated radii (A) for the localization of electronic orbitals in U, eigenenergies (eV) as obtained in an atomic 
DFT calculation, and electron binding energies (eV)'^^. 
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FIG. 6: (Color online), (a) Interatomic screening function for different pairs with the correction scheme proposed as compared 
with previous results, (b) Relative deviation with respect to HF in the internuclear distance corresponding to each energy. 
We compare the results obtained for: (1) the pseudopotentials without semicore (original); (2) with the correction proposed in 
this work (open square symbols); (3) the pseudopotential with semicore states; and (4) the ZBL parametrization. Note that 
considerable deviations in (1) appear for energies of hundreds of eV, and this is considerably improved with (2). ZBL gives 
good description of the high energy regime, but fails for energies smaller than ~0.1 keV. 



